function y = nw_balance(x,k)
 % x is a matrix, n is the time periods, h is the dimension
% the estimated covariance matrix of the columns of x
% output: the S matrix, the covariance matrix of x columns
% only use the balanced rows of x to compute the variance covariance matrix
% does not allow x missing in the middle

% 1. get rid of t's that have nan's
a = isnan(x);
b = find(sum(a,2)>0);
x(b,:) = []; % if any row has any nan's, get rid of the whole row
% 2. Lead-lag covariance
[n,h]=size(x);
if k>0
for i=1:k
    b0=x(1:n-i,:);
    b1=x(i+1:n,:);
    ss_1 = b0'*b1/n; % divide by n, instead of n-i
    ss_2 = b1'*b0/n; % divide by n, instead of n-i
    cov_1{i}=ss_1; % ss_1 = ss_2'
    cov_2{i}=ss_2;
end
end
% 3. Contemporaneous covariance
xx = x'*x/n; % contemporaneous variance-covariance matrix
% 4. Add up lead-lag covariances
if k>0
    i=1;
    y=xx;
    while i<k+1
    y=y+(k-i)/k*cov_1{i}+(k-i)/k*cov_2{i};
    i=i+1;
    end
else
y=xx;
end
end